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ABSTRACT 

Results  from  a  direct  numerical  simulation  of  laminar  flow  over 
a  flat  surface  with  spherical  roughness  elements  using  a  spectral- 
element  method  are  given.  The  numerical  simulation  approximates 
roughness  as  a  cellular  pattern  of  Identical  spheres  protruding  from  a 
smooth  wall.  Periodic  boundary  conditions  on  the  domain’s  horizon¬ 
tal  faeces  simulate  an  inflnite  array  of  roughness  elements  extending 
in  the  streamwise  and  spanwise  directions,  implies  the  parallel-flow 
assumption,  and  results  <n  a  closed  domain.  A  body  force,  designed  to 
yield  the  horizontal  Blasius  velocity  in  the  absence  of  roughness,  sus¬ 
tains  the  flow.  Instabilities  above  a  critical  Reynolds  number  reveal 
negligible  oscillations  in  the  recirculation  regions  behind  each  sphere 
and  in  the  free  stream,  high-amplitude  oscillations  in  the  layer  di¬ 
rectly  above  the  spheres,  and  a  mean  profile  with  an  inflection  point 
near  the  sphere’s  crest.  The  inflection  point  yields  an  unstable  layer 
above  the  roughness  (where  W\y)  <  0)  and  a  stable  region  within 
the  roughness  (where  W{y)  >  0).  Evidently,  the  instability  begins 
when  the  low- momentum  or  wake  region  behind  an  element,  being 
the  region  most  affected  by  disturbances  (purely  numerical  in  this 
case),  goes  unstable  and  moves.  In  incompressible  flow  with  peri¬ 
odic  boundaries,  this  motion  sends  disturbances  to  all  regions  of  the 
domain.  In  the  unstable  layer  just  above  the  inflection  point,  the  dis¬ 
turbances  grow  while  being  carried  downstream  with  a  propagation 
speed  equal  to  the  local  mean  velocity;  they  do  not  grow  amid  the  low 
energy  region  near  the  roughness  patch.  The  most  amplified  distur¬ 
bance  eventually  arrives  at  the  next  roughness  element  downstream, 
perturbing  its  wake  and  inducing  a  global  response  at  a  frequency 
governed  by  the  streamwise  spacing  between  spheres  and  the  mean 
velocity  of  the  most  amplified  layer. 
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1  Introduction 


Generally,  surface  roughness  promotes  transition  in  the  sense 
that  under  otherwise  identical  conditions,  transition  occurs  at  a  lower 
Reynolds  number  on  a  rough  wall  than  on  a  smooth  one  [Dryden,  1953, 
Sedney,  1973].  Evidently,  roughness  elements  give  rise  to  additional 
disturbances  which  add  to  those  already  present  in  the  laminar, 
boundary  layer.  A  sufficiently  rough  surface  may  advance  the  tran¬ 
sition  location  upstream  of  the  smooth-wall  position;  otherwise,  ex¬ 
tended  regions  of  relatively  steady  but  distorted,  laminar  flow  may 
appear  [Sedney,  1973].  Roughness  beginning  downstream  of  the  loca¬ 
tion  where  linear  disturbances  grow  may  enhance  growth  in  Tollmien- 
Schlichting  (TS)  waves  [Corke  &  Others,  1986],  while  disturbances 
from  roughness  beginning  further  upstream  may  “bypass”  the  known 
linear  instability  processes  [Reshotko  &  Leventhal,  1981]. 

1.1  Isolated  roughness 

The  characteristics  of  flow  over  three-dimensional  roughness  ele¬ 
ments  depend  on  the  interaction  between  the  natural  boundary-layer 
vorticity  and  the  obstacle  [Mason  &  Morton,  1987].  Upon  encoun¬ 
tering  a  three-dimensional  surface  protuberance,  the  vortex  lines  may 
concentrate  to  form  a  visible  vortex  core  which  trails  downstream. 
Experiments  at  low  velocity  reveal  a  single  smoke  streak  trailing 
downstream  from  the  element  [Mochizuki,  1961].  At  higher  veloci¬ 
ties,  a  horseshoe  vortex  and  a  pair  of  smoke  spirals  {chimney  vortices 
fixed  in  space  and  due  to  the  fluid  entering  into  the  rear  separation 
pocket)  appear  close  behind  the  sphere  forming  two  parallel  smoke 
filaments  trailing  downstream  parallel  to  the  wall  (called  the  trail¬ 
ing  vortex).  As  the  velocity  increases,  the  trailing  vortex  filaments 
begin  “waving,”  while  growing  over  some  distance  downstream.  Fur¬ 
ther  downstream,  they  cease  to  grow  after  oblique  stretching  by  the 
main  stream’s  shearing  stresses.  At  higher  wind  speed  these  vortices 
shed  at  a  dimensionless  frequency  (non-dimensionalized  as  a  Strouhal 
number,  54  =  fklUk)  ranging  from  0.1  to  0.4  depending  on  Reyn¬ 
olds  number  [Gupta,  1980,  Gregory  &  Walker,  1951].  The  trailing 
vortex  filaments  show  waviness  at  Rck  =  Ukkfu  =  350  and  shedding 
at  Rek  =  700  [Mochizuki,  1961].  (The  undisturbed  velocity  at  the 
crest  of  the  roughness  equals  Uk,  and  the  roughness  height  equals 
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k. )  In  all  cases,  the  dimensionless  instability  frequency,  falls 

above  the  unstable  region  of  the  stability  diagram  for  TS  waves. 

l. 2  Distributed  roughness 

In  contrast  to  the  well  understood  sequence  of  instabilities  lead¬ 
ing  to  transition  around  isolated  roughness,  the  process  for  distributed 
roughness  remmns  unknown  [Morkovin,  1989].  Although  some  simi¬ 
lar  mechanisms  may  exist,  the  ordered  behavior  near  single  elements 
vanishes  when  neighboring  elements  are  present.  We  expect  an  inter¬ 
action  between  the  streamwise  vortices  engendered  by  upstream  elp- 
ments  ana  the  vortex  system  of  those  further  downstream;  however, 
amalgamation  of  the  wakes  from  upstream  elements  should  produce 
a  low-momentum  region  amid  the  roughness  patch.  This  lowers  the 
effective  roughness  Reynolds  number,  Rek,  which  depends  on  the  ve¬ 
locity  at  the  elevation  of  the  roughness  in  the  smooth-wall  boundary 
layer.  Therefore,  with  equal  roughness  Reynolds  numbers,  a  element 
in  a  distributed  roughness  patch  encounters  an  extenuated  main¬ 
stream  compared  to  an  isolated  element.  In  this  case,  closely- packed 
elements  wiU  not  beget  the  strong  horseshoe  and  chimney  vortices 
necessary  to  initiate  the  shedding  process  typically  found  with  iso¬ 
lated  elements  until  reaching  a  much  higher  Reynolds  number. 

2  Numerical  Simulation  of  Roughness 

Several  numerical  studies  investigated  the  mean  flow  and  dis¬ 
turbance  field  around  isolated  three-dimensional  roughness  elements 
[Mason  &  Sykes,  1979,  Mason  &  Morton,  1987,  Tadjfar,  1990).  Tad- 
jfar  employed  a  three-dimensional  triple-deck  with  time-harmonic 
free-stream  disturbances;  and  Mason  et  al.  solved  the  full  three- 
dimensional,  unsteady  Navier-Stokes  equations.  The  work  presented 
here  represents  the  first  computational  investigation  of  distributed 
roughness  using  both  realistic  geometry  and  the  three-dimensional, 
unsteady  Navier-Stokes  equations. 

2.1  Distributed  roughness  pattern  under  study 

The  distributed  roughness  under  investigation  consists  of  an  ar¬ 
ray  of  identical  spheres  protruding  from  a  flat  plate  according  to  pre- 
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Figure  1:  Side  view  of  the  standard  density  domain  showing  a  5X5X5 
mesh  along  the  vertical  symmetry  plane  extending  from  the  wall  to 
the  free  stream. 

vious  experiments  [Kendall,  1981,  Tadjfar  &  Others,  1985].  Whereas 
the  experimental  studies  consisted  of  a  spatial  array  of  a  large  num¬ 
ber  of  spheres,  the  computational  domain  consists  of  the  single  cell 
shown  in  figures  1  and  2. 

2.2  Boundary  conditions 


All  previous  experiments  with  distributed  roughness  were  done 
in  spatially-developing  boundary  layers  or  “open”  flows.  By  applying 
periodic  boundary  conditions  to  the  horizontal  faces  of  a  single  cell 
and  an  impermeable  boundary  condition  along  the  cell’s  top,  we 
simulate  an  infinite  array  of  roughness  elements  extending  upstream, 
downstream,  and  spanwise.  These  boundary  conditions  effectively 
“close”  the  domain  and  result  in  a  parallel  flow. 

In  simple  geometries,  periodic  conditions  allow  us  to  simulate 
spatial  development  by  interpreting  the  time-dependent  data  as  that 
coming  from  a  domain  moving  downstream  with  the  group  velocity 
of  the  disturbances.  This  “moving  box”  approach  is  often  effective 
in  these  geometries  [Singer  &  Others,  1986,  Spalart  &  Yang,  1987, 
Laurien  &  Kleiser,  1989). 

This  artifice  fails  in  complex  geometries,  however.  The  domain 
must  remain  fixed  at  a  given  streamwise  position,  and  temporal  data 
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Figure  2:  Plan  view  of  the  standard  density  domain  showing  a 
5X5X5  mesh  along  a  horizontal  plane  passing  through  the  center 
of  the  spheres. 

can  only  be  used  to  indicate  frequency  amd  phase  content  along  with 
relative  amplitude  between  signals  from  different  locations;  it  cannot 
reveal  anything  about  growth  rates  since  disturbances  cannot  grow  at 
a  fixed  location  in  a  flow  with  steady  mean  quantities.  The  state  with 
steady  mean  flow  is  called  the  time-asymptotic  state.  This  restriction 
does,  however,  provide  a  convenient  check  on  the  numerical  scheme: 
After  achieving  the  time-asymptotic  stage,  a  solution  showing  grow¬ 
ing  fluctuations  indicates  poor  spatial  or  temporal  resolution,  while 
decaying  fluctuations  cannot  occur.  (These  comments  do  not  apply 
to  time- asymptotic  flows  with  imposed  disturbances.  In  these  cases, 
a  time-asymptotic  flow  may  be  “perturbed”  by  introducing  artificial 
disturbances  which  may  grow  or  decay  in  time.) 

Finally,  a  word  about  roughness  and  TS  waves  is  appropriate. 
Although  experiments  reveal  enhanced  growth  in  the  TS  band  of 
frequencies  when  roughness  begins  downstream  of  the  critical  point 
for  TS  wave  growth,  the  periodic  and  free-stream  boundary  condi¬ 
tions  applied  to  the  single  roughness  cell  shown  in  figure  1  preclude 
TS  waves — TS  waves  have  wavelengths  which  typically  equal  ten 
boundary-layer  thicknesses  and  amplitudes  which  decay  asymptoti¬ 
cally  with  distance  from  the  wall.  Unfortunately,  the  finite  height  do¬ 
main,  imposed  free-stream  boundary  conditions,  and  relatively  short 
cell  length  are  incongruous  with  this  type  wave. 
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2.3  Numerical  procedure 


Both  the  incompressible  continuity  equation,  given  by 


duj 

dxj 


=  0, 


and  the  incompressible  Navier-Stokes  equation,  given  by 


dui  dui 


1  dp  ,  d^uj 
pdxi  *  ^  dxjdxj' 


(1) 

(2) 


are  solved  in  the  domain  described  in  sections  2.1  and  2.2  using 
a  spectral-element  code.  The  spectral-element  method  employs  a 
three-step,  time-splitting  scheme  where  the  non-linear,  pressure,  and 
viscous  terms  in  the  incompressible  Navier-Stokes  equations  are  writ¬ 
ten  in  separate  fractional  steps  [Orszag  &  KeUs,  1980].  By  introduc¬ 
ing  intermediate  velocities,  these  steps  yield  consecutive  solutions 
based  first  on  the  non-linear,  then  on  the  pressure,  and  finally  on  the 
viscous  terms. 


2.4  Input  parameters 


When  the  incompressible  Navier-Stokes  equations  (see  equa¬ 
tion  2)  are  written  in  non-dimensional  form,  three  dimensionless 
quantities  appear:  the  Reynolds  number.  Re  =  ULIv;  the  StrouhaJ 
number,  St  =  fL/U’,  and  the  Froude  number,  F,  =  U/y/gL.  With 
periodic  boundary  conditions,  a  body  force  generates  the  velocity, 
and  a  relationship  exists  between  the  Reynolds  and  Froude  num¬ 
bers.  Furthermore,  in  a  flow  without  external  forcing,  the  Strouhal 
number  governs  unsteadyness  due  to  naturally  arising  instabilities. 
Hence,  the  Reynolds  number  or,  for  a  given  geometry,  the  fluid  vis¬ 
cosity  and  the  body  force,  completely  determines  the  flow. 

Here  we  employ  a  special  body  force  which,  when  twice  differ¬ 
entiated  with  respect  to  y,  yields  the  Blasius  profile  in  the  absence 
of  roughness  elements.  It  is  derived  by  considering  the  streamwise 
momentum  equation  in  a  boundary  layer  on  a  smooth,  flat  surface 
with  no  imposed  streamwise  pressure  gradient.  Equation  2  in  this 
case  reduces  to 


du  du  d^u 


(3) 
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By  interpreting  the  right-hand  side  of  equation  3  as  a  body  force,  we 


get 

du  du  ... 

“ai  +  "ai  = 

(4) 

with 

(5) 

while  the  subscript  B  refers  to  the  Blasius  value.  By  twice  differenti¬ 
ating  the  Blasius  velocity  profile  according  to  equation  5,  we  obtain 
a  body  force  which  yields  the  Blasius  profile  on  a  smooth,  flat  surface 
with  no  imposed  pressure  gradient.  When  this  body  force  is  applied 
to  a  domain  with  surface  roughness  elements,  the  resulting  velocity 
Held  represents  the  influence  of  the  roughness  elements  on  the  Blasius 
boundary  layer. 

2.5  Time  asymptotic  state 


From  a  given  initial  state,  a  certain  time  elapses  while  mean- 
flow  quantities  develop  towards  a  time-asymptotic  state.  The  time 
required  to  reach  this  state  depends  on  the  starting  solution  as  well 
as  the  longer  of  the  convection  and  diffusion  times. 

The  convection  time  to  reach  the  asymptotic  state  depends  on  the 
time  for  a  particle  to  pass  across  the  domain,  given  by  =  LxjUoo, 
where  Uoo  equals  the  streamwise  free-stream  velocity,  and  Lx  equals 
a  streamwise  length  scale. 

The  diffusion  time  to  reach  the  time- asymptotic  state  depends 
on  momentum  diffusion  across  the  layer.  For  a  fixed  boundary-layer 
thickness  and  viscosity,  momentum  diffusion  across  a  boundary  layer 
of  thickness  S  requires  time. 

The  ratio  of  diffusion  to  convection  times  equals 


Td 

Tc 


fi^lv 

LxlU^ 


,  6  .UooS  ,  6  -p 


(6) 


Typically,  the  Reynolds  number  equals  several  hundred,  and  the 
ratio  of  boundary-layer  thickness  to  streamwise  length  is  of  order 
one;  hence,  diffusion  requires  two  orders-of-magiutude  more  time 
than  convection.  For  typical  flows  of  interest,  diffusion  requires 
10^  time  units.  Since  the  first-order  time-accurate  spectral-element 
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Figure  3:  Side  view  of  the  high-deusity  domain  showing  a  5X5X5 
mesh  along  the  vertical  symmetry  plane  extending  from  the  wall  to 
the  free  stream. 

method  requires  small  time  steps — typically  At  <  0.01 — for  good  ac¬ 
curacy,  approximately  10®  time  steps  elapse  before  reaching  a  time- 
asymptotic  flow. 

2.6  Cases  studied 


Some  previous  experiments  using  distributed  roughness  elements 
[Gartshore  &  de  Croos,  1979,  Lee  &  Soliman,  1977,  Kendall,  1981] 
found  that  roughness  “density,”  oi  spacing  between  elements,  con¬ 
trols  certain  mechanisms  responsible  for  transition.  Based  on  this 
idea,  I  selected  three  different  roughness  densities  for  study:  stan¬ 
dard,  low,  and  high.  The  roughness  density  was  varied  by  stretching 
or  compressing  the  domain  in  the  streamwise  direction  while  holding 
the  spanwise  dimension,  the  number  of  finite  elements,  aind  the  to¬ 
tal  number  of  grid  points  constant.  Consequently,  the  high-density 
domain  was  shorter  (see  figure  3),  and  the  low-density  domain  was 
longer  (see  figure  4)  than  the  standard  domain.  The  standard-density 
domain  consists  of  a  square  with  streamwise  and  spanwise  lengths  of 
8.4;  the  high-density  domain  consists  of  a  rectangle  with  a  streamwise 
length  of  6.0  and  spanwise  length  8.4;  and  the  low-density  domain 
has  a  streamwise  length  of  12.0  and  a  spanwise  length  of  8.4.  All 
cases  used  Rck  =  175,  k/6*  =  0.72,  and  the  Blasius  based  body  force 


Figure  4:  Side  view  of  the  low-density  domain  showing  a  5X5X5 
mesh  along  the  vertical  symmetry  plane  extending  from  the  wall  to 
the  free  stream. 

(see  section  2.4)  to  drive  the  flow.  For  each  case,  the  one-cell  grid 
consisted  of  144  elements  and  a  5Jf5X5  mesh  within  each  clement. 
The  standard-density  case  was  also  run  with  a  finer,  7X7X7,  mesh. 

The  calculations  performed  on  a  single-cell  domain  correspond  to 
the  assumption  of  “strong  periodicity”  [Karniadakis  &  Others,  1988]. 
Unfortunately,  we  have  no  a  priori  knowledge  whether  or  not  a  given 
flow  will  behave  differently  when  the  number  of  cells  are  changed.  To 
gain  insight  into  the  effect  of  the  streamwise,  periodic  boundary  con¬ 
ditions  on  the  nature  of  the  disturbance  wavelength  and  frequency,  a 
two-cell  high-density  case  and  a  two-ceU  low-density  case  were  run. 
These  two-cell  cases  contained  two  144  -element  grids  aligned  in  the 
streamwise  direction,  used  the  5X5X5  mesh,  and  were  twice  as  long 
as  their  single-cell  counterparts.  A  change  in  results  between  the 
one-  and  two-ceU  cases  will  indicate  the  presence  of  a  sub-harmonic 
disturbance. 

3  Steady  State  Results 


Since  we  compute  the  unsteady  Navier- Stokes  equations  until 
a  time-asymptotic  state  is  reached,  all  of  the  results  are  unsteady 
above  a  critical  Reynolds  number.  To  distinguish  steady  results  from 
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unsteady  results  (given  in  the  next  section)  the  term  “steady  state” 
is  used  loosely  to  refer  to  quantities  which  do  not  vary  significantly 
with  time.  These  Include  mean  flow  boundary-layer  quantities  and 
instantaneous  velocity  and  vorticity  contours. 

3.1  Mean  flow  results 

Results  for  the  three  different  roughness  densities  using  a  5X5X5 
mesh  with  one  cell  are  given  in  table  1.  Also  appearing  in  table  1 
are  results  for  the  standard-density  case  using  a  7X7X7  mesh.  The 
corresponding  smooth-wall,  Blasius  results  are  given  for  comparison. 
Some  of  the  quantities  appearing  in  the  table  need  to  be  defined  be¬ 
fore  discusing  the  results.  The  ratio  of  rough-wall  to  smooth-wad 
friction  coefficients  equals  The  maximum  velocity  divergence, 

equals  (V-V)m,  and  the  velocity  divergence  integrated  over  the  entire 
domain,  equals  (V  •  V)a — both  indicate  the  quality  of  the  numerical 
solution,  since  they  equal  zero  in  incompressible  flow.  Finally,  the 
displacement  and  momentum  thicknesses  appear  here  in  dimension¬ 
less  form. 

According  to  table  1,  an  improvement  in  grid  resolution  from 
a  5X5X5  to  a  7X7X7  mesh  per  element  yields  a  slight  increase  in 
drag,  evidence  that  the  course  mesh  fails  to  capture  the  smaller  scales 
required  to  resolve  the  flow.  The  results  also  show  an  increase  in  6’ 
and  a  decrease  in  0  compared  to  the  smooth-wall  (Blasius)  values 
and  consistent  with  inflected  or  adverse  pressure-gradient  profiles. 


density 

mesh 

B 

0 

H 

(V-V^)„ 

stand. 

5X5X5 

1.85 

.508 

3.64 

0.0200 

-0.0176 

2.62 

stand. 

7X7X7 

1.88 

.530 

3.57 

— 

— 

2.64 

high 

5X5X5 

1.80 

.473 

3.80 

1 9 

2.69 

low 

5X5X5 

1.84 

.529 

3.48 

■111  wW 

2.49 

Bias. 

theory 

1.72 

.664 

2.59 

0.0 

0.0 

1.0 

Table  1:  Mean  flow  results  for  the  various  roughness  densities 
and  the  corresponding  smooth-wall,  Blasius  data 


Compared  to  the  standard  domain,  the  high-density  case  effects 
an  increase  in  and  H  and  a  decrease  in  both  6*  and  0,  the 
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Figure  5;  Side  view  of  the  streamwise  velocity  contours  on  a 
standard-density  domain  with  a  SXSJf  5mesh  along  the  vertical  sym¬ 
metry  plane  extending  from  the  wall  to  the  free  stream. 

opposite  occurs  in  the  low-density  case. 


3.2  Velocity  contours  standard  density  domain 

The  standard-density  domain  velocity  contours  are  given  in  fig¬ 
ures  5  through  10;  vorticity  contours  are  given  in  figure  11.  Figures  5 
and  6  show  the  streamwise  and  vertical  velocity  on  a  plane  through 
the  center  of  the  middle  sphere.  A  small  reverse  flow  appears  in  the 
region  behind  the  sphere.  The  vertical  velocity  (figure  6)  equals  zero 
above  two  roughness  heights  while  remaining  within  4.4  percent  of 
the  free-stream  velocity  near  the  elements.  The  small  upward  veloc¬ 
ity  behind  the  element  remains  confined  to  the  region  between  the 
centerline  and  crest.  At  this  Reynolds  number,  the  sign  and  location 
accede  with  the  chimney  vortex  of  an  isolated  roughness  element — 
a  vortex  rotating  counter  to  and  above  the  horseshoe  vortex.  In 
this  case,  however,  the  neighboring  roughness  elements  beget  a  low- 
momentum  region  and  prevent  the  formation  of  robust  horseshoe  and 
chimney  vortices.  Instead,  a  flow  dominated  by  kinematics  results: 
continuity  in  concert  with  the  solid-wall  boundary  conditions.  The 
upward  velocity  behind  the  sphere  depends  on  continuity:  the  fluid 
passing  around  the  sides  of  the  element  converges  in  the  wake  and 
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Figure  6:  Side  view  of  the  vertical  velocity  contours  on  a  standard- 
density  domain  with  a  5X5X5  mesh  along  the  vertical  symmetry 
plane  extending  from  the  wall  to  the  free  stream. 

gets  forced  upward  (the  wall  prevents  it  from  going  down). 

Reverse-flow  regions  indicate  the  importance  of  inertia  since  no 
reverse-flow  regions  exist  in  a  viscous  dominated  flow.  Thus,  the 
upward  and  upstream  velocity  in  the  wake  indicates  that  a  chimney 
vortex,  albeit  a  weak  one,  forms.  Does  the  horseshoe  vortex  ap¬ 
pear?  With  an  isolated  sphere,  the  horseshoe  vortex  begins  in  the 
stagnation  region  and  extends  into  the  wake.  In  distributed  rough¬ 
ness,  however,  other  roughness  elements  disturb  the  velocity  in  the 
wake;  consequently,  the  stagnation  region  becomes  the  best  place  to 
evince  the  horseshoe  vortex.  Unfortunately,  only  an  extremely  weak 
horseshoe  vortex  exists.  A  true  horseshoe  vortex  would  appear  in  the 
stagnation  region  as  a  downward  velocity  near  the  element’s  center- 
line,  an  upstream  velocity  near  the  wall,  and  an  upward  velocity  out 
in  front  of  the  element.  As  shown  in  figure  6,  the  downward  velocity 
occurs  near  the  sphere;  however,  neither  the  upward  nor  upstream 
velocities  are  evident. 

The  same  velocity  components  on  a  horizontal  plane  passing 
through  the  center  of  the  spheres  are  shown  in  figures  7  and  8.  A  high 
streamwise  velocity  “channel”  appears  between  the  rows  of  spheres 
as  a  consequence  of  forcing  the  flow  along  the  geometric  symmetry 
plane,  and  a  reverse-flow  region  appears  downstream  of  each  ele¬ 
ment.  The  vertical  velocity  contours  show  downward  flow  ahead  and 
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Figure  7:  Plan  view  of  the  stream  wise  velocity  contours  on  a 
standard-density  domain  with  a  5X5X5  mesh  along  a  horizontal 
plane  passing  through  the  center  of  the  spheres. 


Figure  8:  Plan  view  of  the  vertical  velocity  contours  on  a  standard- 
density  domain  with  a  5X5X5  mesh  along  a  horizontal  plane  passing 
through  the  center  of  the  spheres. 
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Figure  9:  Cross-sectional  view  of  the  vertical  velocity  contours  on  a 
standard-density  domain  with  a  5X5X5  mesh  along  a  vertical  plane 
passing  through  the  domain’s  inlet  extending  from  the  wall  to  the 
free  stream. 

upward  flow  behind  each  element,  consistent  with  the  continuity  ar¬ 
gument  previously  outlined. 

Figures  9  and  10  present  the  vertical  and  spanwise  velocity  on  a 
cross-sectional  plane  along  the  inlet  of  the  domain.  This  plane  cuts 
through  the  center  of  the  corner  spheres  and  lies  just  downstream 
of  the  center  sphere  in  an  imaginary  cell  upstream  of  the  computa¬ 
tional  cell.  Hence,  this  plane  shows  the  central  sphere’s  wake  struc¬ 
ture.  The  spanwise  velocity  always  points  toward  the  centerline, 
clearly  showing  how  the  fluid  fills  in  the  region  behind  the  center 
sphere  while  being  pushed  inward  by  the  corner  spheres — both  ef¬ 
fects  act  in  the  same  direction,  yielding  a  strong  spanwise  velocity. 
The  vertical  velocity  pattern  illicits  more  interest.  The  central  re¬ 
gion  shows  upward  velocity  consistent  with  the  continuity  constraint 
discussed  previously.  A  larger  upward  velocity  region  develops  above 
each  sphere — a  manifestation  of  the  inertia  of  the  fluid  passing  over 
the  obstacle.  (If  the  flow  were  free  of  inertia,  all  velocity  on  this  plane 
would  be  in  the  streamwise  direction.)  The  two  downward  velocity 
regions  reflect  wake  filling  from  above  the  corner  spheres. 

Streamwise  vorticity  on  a  cross-sectional  plane  at  the  inlet  is 
shown  in  figure  11.  These  vorticity  regions  seem  to  match  those  pro¬ 
duced  by  the  horseshoe  and  chimney  vortices  behind  isolated  rough- 
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Figure  10:  Cross-sectional  view  of  the  spanwise  velocity  contours 
on  a  standard-density  domain  with  a  5X5X5  mesh  along  a  vertical 
plane  passing  through  the  domsdn^s  inlet  extending  from  the  wall  to 
the  free  stream. 
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Figure  11:  Cross-sectional  view  of  the  streamwise  vorticity  contours 
on  a  standard-density  domain  with  a  5X5X5  mesh  along  a  vertical 
plane  passing  through  the  domain’s  inlet  extending  from  the  wall  to 
the  free  stream. 
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Figure  12:  Side  view  of  the  vertical  velocity  contours  on  a  high- 
density  domain  with  a  5X5X5  mesh  along  the  vertical  symmetry 
plane  extending  from  the  wall  to  the  free  stream. 

ness  elements.  Congruent  with  the  vorticity  system  downstream  of 
the  center  sphere  without  neighbors,  on  the  left  we  see  positive  vor¬ 
ticity  below  and  negative  vorticity  above  the  midpoint;  the  signs 
reverse  on  the  right.  We  could  easily  link  the  trailing  vortices  from 
the  sphere  upstream  and  the  vorticity  pattern  on  the  inlet  plane. 
However,  if  these  regions  actually  contain  vortices,  circular  velocities 
should  appear  in  the  spanwise  and  vertical  contour  plots.  We  see  no 
evidence  of  circular  motion. 


3.3  Velocity  contours  high  and  low  density  domains 

The  vertical  velocity  contours  on  the  symmetry  plane  (see  fig¬ 
ures  12  and  13  )  show  that  the  dense-packing  maximum  velocity 
decreases  and  the  sparse-packing  value  increases  compared  to  the 
standard  domain.  Evidently  the  sparse  array  allows  an  increase  in 
momentum  between  the  wake  of  one  sphere  and  the  stagnation  region 
of  the  next.  Further  confirmation  appears  in  the  streamwise  velocity 
contours  on  the  horizontal  plane  passing  through  the  sphere  cen¬ 
ters  (not  shown).  The  maximum  streamwise  velocity  equals  0.0875 
in  the  standard  domain,  0.054  in  the  dense  domain,  and  0.115  in 
the  sparse  domain,  a  difference  of  a  factor  of  two  between  the  dense 


16 


-40zrs 

•00325 

•00173 

-00125 

-00073 


00075 

00125 

00175 

00325 

00273 

00325 

Umk 


Figure  13:  Side  view  of  the  vertical  velocity  contours  on  a  low-density 
domain  with  a  5X5X5  mesh  along  the  vertical  symmetry  plane  ex¬ 
tending  from  the  wall  to  the  free  stream. 

and  sparse  cases.  Clearly  the  dense  array  creates  a  large  region  of 
low-momentum  fluid  amid  the  spheres.  An  even  more  striking  ef¬ 
fect  is  evident  in  the  span  wise  and  vertical  velocity  contours  on  the 
inlet  cross-sectional  plane  (figures  14  and  15).  Figure  14  shows  the 
spanwise  velocity  contours  for  the  dense  array.  The  usual  conver¬ 
gent  velocity  regions  are  shown;  however,  an  outward  velocity  below 
this  region  exists,  a  downward  velocity  between  these  regions  ap¬ 
pears,  and  an  upward  velocity  outside  the  downward  velocity  region 
develops — the  horseshoe  vortex.  There  is  no  evidence  of  the  chim¬ 
ney  vortex.  Somehow  the  close  proximity  of  the  corner  spheres  to 
the  wake  of  the  center  sphere  generated  the  horseshoe  vortex. 


4  Time  Dependent  Results 


During  each  run,  we  store  streamwise,  spanwise,  and  vertical 
velocity  signals  from  five  numerical  probes.  Probe  locations  are  given 
in  figures  16  and  17.  Probe  number  48,  located  behind  the  center 
sphere,  lies  within  one  sphere  elevation  from  the  wall;  probe  numbers 
62,  78,  and  88  lie  approximately  one  and  one-half  sphere  diameters 
from  the  wall;  and  probe  number  120  lies  two  and  one-half  diameters 
above  the  wall.  Probe  numbers  62  and  78  occupy  symmetric  positions 
behind  each  of  the  two  forward  corner  spheres.  Similarly,  probes  78 
and  88  occupy  equivalent  positions  behind  the  center  sphere  and  one 
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Figure  14:  Cross-sectional  view  of  the  spajiwise  velocity  contours  on 
a  high-density  domain  with  a  5X5X5  mesh  along  a  vertical  plane 
passing  through  the  domain’s  inlet  extending  from  the  wall  to  the 
free  stream. 
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Figure  15:  Cross-sectional  view  of  the  vertical  velocity  contours  on 
a  high-density  domain  with  a  5X5X5  mesh  along  a  vertical  plane 
passing  through  the  domain’s  inlet  extending  from  the  wsJl  to  the 
free  stream. 
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Figure  16:  Side  view  of  the  probe  locations. 


Figure  17:  Plan  view  of  the  probe  locations. 


19 


Figure  18:  Vertical  (top  figure),  streamwise  (middle  figure),  and 
spanwise  (bottom  figure)  velocity  versus  time  for  the  five  probes  in 
the  standard-density  domain. 


of  the  forward  corner  spheres.  Thus,  differences  in  the  signals  from 
each  of  these  probes  may  reveal  flow  asymmetries  and  phase  shifts 
in  the  disturbance  waves. 


4.1  Standard  density  domain 


Using  the  standard  domain,  mesh  dependent  solutions  revealed 
both  steady  and  unsteady  time-asymptotic  states.  Unsteady  flows 
using  a  5X5X5  mesh  decayed  with  a  7X7X7  mesh. 

Results  using  a  5X5X5  mesh  are  given  in  figure  18.  Although 
not  shown,  the  total  kinetic  energy  variation  (an  indication  of  flow 
development  toward  the  time-asymptotic  state)  was  less  than  0.1 
percent  over  the  time  contained  in  the  graphs.  This  corresponds  to 
a  dimensionless  convection  time,  Tc  =  ^tUoolLx,  of  66,  or  a  dimen¬ 
sionless  diffusion  time,  equal  to  0.51.  These  numbers 

indicate  that  a  fluid  element  in  the  free  stream  passed  across  the 
domain  66  times,  but  momentum  only  diffused  across  approximately 
0.7A:  during  the  measuring  interval.  The  small  change  in  total  kinetic 
energy,  however,  indicates  that  the  flow  may  appropriately  classified 
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Figure  19:  Vertical  velocity  amplitude  for  each  of  the  five  probes 
versus  Strouhal  number  in  the  standard-  (top  figure),  high-  (middle 
figure),  and  low-density  (bottom  figure)  domains. 

as  time-asymptotic.  No  oscillations  developed  in  any  of  the  velocity 
components  from  probe  numbers  48  and  120  until  the  r.m.s.  value  of 
probe  78’s  streamwise  velocity  reached  0.5  percent  of  the  free-stream 
velocity.  At  this  point  probe  120  began  a  low-level  oscillation.  Ver¬ 
tical,  streamwise,  and  spanwise  oscillations  with  a  180°  phase  shift 
occurred  in  probes  78  and  88;  a  small  phase  shift  exists  between 
probes  62  and  88.  These  oscillations  produce  a  strong  single  spike  at 
Stk  =  kffUk  =  0.33  in  the  standard  domain’s  frequency  spectrum 
(see  figure  19).  This  corresponds  to  the  “passing”  frequency,  UkILx, 
for  the  time  a  parcel  of  fluid  at  elevation  y  =  k  requires  to  cross 
between  spheres  spaced  at  a  distance  Lx  with  a  Blasius  velocity  Uk- 

4.2  High  and  low  density  domains 

If  the  roughness  density  controls  the  characteristic  frequency,  a 
lower  oscillation  frequency  would  exist  in  a  low- density  domain,  and 
a  higher  one  in  a  high-density  domain.  The  high-density  domain’s 
streamwise  length  equals  6.0;  the  low-density  domain’s  streamwise 
length  equals  12.0.  These  lengths  correspond  to  passing  frequency 
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Figure  20:  Vertical  (top  figure),  streamwise  (middle  figure),  and 
spanwise  (bottom  figure)  velocity  versus  time  for  the  five  probes  in 
the  high-density  domain. 

based  Strouhal  numbers  of  0.47  and  0.23  for  the  high-  and  low-density 
domains  respectively.  The  time- dependent  data  from  the  five  probes 
are  given  in  figures  20  and  21  for  these  cases.  Although  not  obvious 
from  these  figures,  the  dominant  oscillation  frequency,  according  to 
the  frequency  spectra  shown  in  figure  19,  does  agree  with  the  pre¬ 
dicted  Strouhal  numbers  given  above. 

5  Analysis 

5.1  Mean  flow  variations  with  roughness  density 

At  first  glance  it  would  appear  that  the  displacement  and  mo¬ 
mentum  thicknesses  would  increase  as  the  roughness  elements  be¬ 
come  more  densely  packed.  However,  in  the  limit  of  infinite  roughness 
density,  these  quantities  would  equal  their  Blasius  values  since  the 
roughness  effectively  becomes  a  flat  plate.  A  plot  of  6*  and  0  versus 
roughness  density  would  show  Bl<isius  values  at  both  the  infinite- 
density  and  zero-density  limits.  In  between,  the  displacement  thick¬ 
ness  increases  and  the  momentum  thickness  decreases  compared  to 
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Figure  21:  Vertical  (top  figure),  streamwise  (middle  figure),  and 
spanwise  (bottom  figure)  velocity  versus  time  for  the  five  probes  in 
the  low-density  domain. 

Blasius  values,  consistent  with  inflected  or  adverse  pressure  gradient 
profiles. 

Another  parameter  which  changes  with  roughness  density  is  the 
friction  coefl^cient.  Schlichting  investigated  how  the  flow  resistance 
depends  on  roughness  density  [Schlichting,  1936].  He  defined  rough¬ 
ness  density  as 

^  _  projected  area  of  roughness  on  plate 
^  plate  area 

and  found  that  the  maximum  resistance  does  not  occur  at  the  max¬ 
imum  roughness  density  but  at  a  considerably  lower  value: 

{Fr)max  =  0.4.  (8) 

The  three  roughness  densities  used  in  this  study  correspond  to  /v  = 
0.20  in  the  standard  domain,  Fr  =  0.28  in  the  high-density  domain, 
and  Fr  =  0.14  in  the  low-density  domain.  Schlichting  defines  a  re¬ 
sistance  coefficient  in  terms  of  the  additional  drag  due  solely  to  the 
roughness  elements  as  Wr  and  the  corresponding  friction  coefficient 
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due  to  the  roughness  elements  as 


f  trough  only  — 


Wr 


\PVIF,- 


(9) 


Defining  the  ratio  of  rough-  to  smooth-wall  friction  factors  as  x  gives 

//->  ■>  s)tmootk{X  ~  f)  /inx 

f  jTough  only  —  p  ■  (1") 

Using  the  given  values  of  iv  and  the  results  from  table  1  for  x  gives 
the  results  shown  in  table  2.  Also  included  are  Schlichting’s  experi¬ 
mental  results. 


Case 

mm 

Uk/Uoo 

standard-density 

iBI 

0.43 

high-density 

0.43 

low-density 

iSI 

0.43 

Schlichting 

0.126 

0.704 

Schlichting 

0.349 

0.678 

0.017 

Schlichting 

0.907 

0.826 

0.007 

Schlichting 

0.126 

0.570 

0.011 

Table  2:  Roughness  friction  coefficient  versus  roughness  density, 
comparison  between  numerical  results  and  experiments. 


Although  the  smallest  total  friction  factor  occurs  in  the  low- 
density  domain,  we  see  that  the  friction  factor  due  to  the  spheres 
alone  is  greatest  in  this  case.  The  low-density  domain  permits  more 
acceleration  between  spheres  resulting  in  a  faster  approaching  stream; 
consequently,  each  element  produces  greater  resistance. 

The  numerical  cases  studied  laminar  flow  and  yielded  rough¬ 
ness  resistance  coefficients  higher  than  Schlichting’s  turbulent-flow 
results.  This  comports  with  expected  behavior:  an  increase  in  fric¬ 
tion  coefficient  with  decreasing  Reynolds  number. 

5.2  Rough  to  smooth  wall  comparison 


The  previous  discussion  concerns  the  effect  of  roughness  den¬ 
sity  on  the  rough-wall  bound  ary- layer  integral  quantities;  this  sec¬ 
tion  relates  the  roughness  quantities  to  the  corresponding  Blasius  or 
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smooth-wall  values.  All  cases  reveal  a  larger  dimensiooless  displace¬ 
ment  thickness  than  the  Blasius  value  1.72,  a  dimensionless  momen¬ 
tum  thickness  smaller  than  the  Blasius  value  0.664,  and  a  shape 
factor  larger  than  the  Blasius  value  2.59.  These  results  reflect  the 
nature  of  the  roughness  boundary  layer:  a  low-momentum  region 
surrounding  the  roughness  elements  and  a  high-momentum  region 
just  above  the  roughness.  Evidently,  the  spheres  push  fluid  toward 
the  free  stream,  leaving  the  region  within  one  diameter  from  the  wall 
with  low  momentum.  Near  the  wall,  the  term  (1  —  U/Uoo)  is  large, 
yielding  a  larger  6*.  Also  near  the  wall,  UjUoo  is  small,  yielding  a 
lower  $. 

The  difference  between  the  rough-  and  smooth-wall  displacement 
thicknesses  indicates  the  extent  of  the  velocity  profile  displacement 
due  to  roughness.  For  the  smooth  wall  we  can  estimate  the  dis¬ 
placement  thickness  as  as  ^/3  =  3.67;  a  typical  rough- wall  mea¬ 
sured  value  gives  6*  as  3.90.  Therefore,  the  additional  displace¬ 
ment  due  to  roughness  normalized  by  the  roughness  height  equals 
A6*/k  =  (3.90  — 3.67)/2.8  =  0.082.  Comparing  this  measured  change 
in  displacement  thickness  to  the  volumetric  average  thickness  of  the 
roughness  itself, 


Ah 

k 


kL^L, 


0.11, 


(11) 


shows  that,  roughness  displaces  the  Blasius  layer  an  additional  three 
quarters  of  the  average  roughness  height,  AS*/ Ah  as  0.75. 


5.3  Streamwise  vorticity 


The  experiments  with  isolated  roughness  show  the  importance  of 
trailing  vortices  on  boundary-layer  instabilities.  With  this  in  mind,  it 
is  of  interest  to  examine  the  distributed  roughness  flow  for  streamwise 
vortices.  With  the  aid  of  figures  9  and  10,  the  spanwise  and  verti- 
c«d  velocity  contours  show  no  evidence  of  circular  flow.  Instead,  the 
vorticity  patterns  depend  on  the  gradients  of  mean  velocity.  Mathe¬ 
matically,  the  streamwise  vorticity  component  equals 

Qx  = -{dxu/dy  -  dv/dz).  (12) 

In  the  results  with  spherical  roughness,  the  variation  of  spanwise 
velocity  with  vertical  elevation,  dw/dy,  dominates  the  other  term. 
Referring  to  figure  10,  we  see  that  w  grows  in  magnitude  from  zero 
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at  the  wall  to  a  maximum  at  the  sphere’s  midplane  and  back  to  zero 
above  two  diameters.  On  the  right  side  of  the  plane,  w  is  positive; 
hence,  dw/dy  >  0  near  the  wall,  and  the  vorticity  is  negative.  Above 
the  midplane,  dwfdy  <  0,  and  the  vorticity  is  positive.  The  opposite 
signs  correspond  to  the  left  half  of  the  plane.  This  pattern  exactly 
matches  the  streamwise  vorticity  regions  of  figure  11.  Evidently,  the 
presence  of  other  spheres  and  the  relatively  low  Reynolds  number 
used  in  the  computations  prevent  the  strong  streamwise  vortices  ob¬ 
served  in  experiments  with  isolated  roughness. 

5.4  Oscillation  frequency  and  phase 

The  instability  scenario  begins  when  the  local  recirculation  re¬ 
gions  behind  each  sphere  become  unsteady  above  a  certain  critical 
Reynolds  number  (still  unknown).  This  unsteady  motion  affects  the 
mean  flow  by  introducing  pressure  and  velocity  disturbances  which, 
due  to  the  periodic  boundary  conditions  and  incompressibility  of  the 
fluid,  instantly  fill  the  entire  domain. 

Once  introduced,  these  disturbances  grow  or  decay  while  prop¬ 
agating  downstream  with  the  local  convection  velocity  of  the  mean 
flow.  Eventually,  those  disturbances  receiving  the  most  energy  in 
unstable  modes  dominate  the  flow.  It  appears  that  growth  does  not 
occur  in  the  recirculation  regions  behind  each  sphere,  since  there 
is  little  incentive  for  amplification  this  close  to  the  wall.  Instead, 
the  major  growth  occurs  in  the  destabilized  region  just  above  the 
inflection  point. 

This  comports  with  the  stability  characteristics  governed  by  U"{y). 
Inflected  profiles  usually  have  a  higher  amplification  rate  for  the  same 
frequency  and  Reynolds  number  than  a  profile  with  purely  negative 
U"(y)  [Obremski  &  Others,  1969].  Although  inflected  profiles  are 
less  stable,  there  is  a  stabilizing  effect  in  the  region  between  the  in¬ 
flection  point  and  wall  where  U"{y)  >  0.  An  examination  of  the 
streamwise  momentum  equation  reveals  that  U'\y)  acts  like  a  force 
per  unit  mass  on  a  parcel  of  fluid.  When  positive,  as  within  the 
region  between  the  wall  and  the  inflection  point,  it  accelerates  or 
stabilizes  a  fluid  parcel;  outside  of  the  inflection  point  it  opposes 
or  destabilizes  the  motion.  This  explains  why  oscillations  at  probe 
48,  located  below  the  inflection  point,  are  small;  why  oscillations  at 
probes  62,  78,  and  88,  located  just  above  the  inflection  point,  are 
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large;  and  why  oscillations  at  probe  120,  located  in  a  region  where 
U"{y)  is  small,  are  small.  What  explains  the  oscillation  frequency? 

Apparently,  the  most  amplified  disturbances  occur  in  the  layer 
just  above  the  spheres,  presumably  carried  downstream  with  a  ve¬ 
locity  of  the  local  mean  flow.  Upon  reaching  the  exit  boundary,  the 
disturbances  reappear  along  the  inlet  boundary  then  resume  their 
journey  downstream.  Once  again,  they  encounter  the  exit  bound¬ 
ary,  and  the  pattern  repeats  ad  infinitum.  A  disturbance  interacts 
with  the  low-momentum  recirculation  zones  behind  a  sphere  each 
time  it  passes  overhead.  (It  also  interacts  with  higher  momentum 
regions,  but  these  regions  do  not  respond  significantly.)  If  the  flow 
is  unstable,  the  interaction  increases  the  amplitude  of  the  distur¬ 
bance  and  causes  the  recirculation  zone  to  move  more  vigorously. 
This  generates  a  periodic  motion  with  period  equal  to  the  time  for 
a  fluid  element  to  pass  between  spheres  at  the  elevation  of  the  most 
amplified  disturbance — in  this  case  the  period  equals  LxjUk-  An  in¬ 
dividual  recirculation  region  interprets  this  as  a  periodic  forcing  or 
“buffeting.”  If  this  buffeting  frequency  falls  near  the  body’s  natural 
frequency,  a  response  at  the  buffeting  frequency  occurs.  This  is  re¬ 
ferred  to  as  the  “lock-in”  mode  [Karniadakis  &  Triantafyllou,  1989]. 
When  the  buffeting  frequency  differs  greatly  from  the  natural  fre¬ 
quency,  the  flow  responds  at  the  isolated-body  natural  frequency, 
referred  to  as  “non-lock-in.”  Evidently,  the  three  roughness  density 
cases  under  consideration — standard,  high,  and  low — all  satisfy  the 
lock-in  mode  criteria.  In  each  case,  the  oscillation  frequency  equals 
the  forcing  frequency. 

Presumably,  a  higher  or  lower  density  case  would  show  the  non¬ 
lock-in  mode.  In  this  case,  the  oscillation  frequency  would  equal 
the  isolated  element  natural  frequency  and  would  be  independent 
of  periodic  boundary  conditions  or  element  spacing.  We  may  also 
expect  a  shear-layer  based  instability  instead  of  one  based  on  vortex 
shedding. 

According  to  Karniadikis  and  Triantafyllou,  1989,  the  distur¬ 
bance  or  forcing  amplitude  also  influences  the  response.  They  present 
a  plot  of  forcing  amplitude  versus  forcing  frequency,  with  unsteady 
regions  labeled  “lock-in”  and  “non-lock-in.”  The  lock-in  region  oc¬ 
curs  within  an  open-end-up  parabolic  curve  centered  at  the  natural 
frequency  of  the  isolated  body.  The  non-lock-in  region  surrounds  the 
lock-in  zone.  Furthermore,  the  unsteady  region — whether  lock-in  or 
non-lock-in — only  occurs  above  a  threshold  amplitude.  Thus,  forcing 
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at  the  natural  frequency  with  an  amplitude  just  above  the  threshold 
will  produce  a  lock-in  response;  forcing  with  the  same  amplitude  at 
a  slightly  different  frequency  will  produce  a  non-lock-in  response.  As 
the  forcing  amplitude  increases,  the  lock-in  mode  region  grows;  pre¬ 
sumably,  at  infinite  amplitude  any  forcing  frequency  will  yield  the 
lock-in  behavior. 

In  the  numerical  calculations  shown  here,  the  only  forcing  is  due 
to  numerical  approximations,  and  it  seems  reasonable  to  assume  that 
the  forcing  level  is  greater  with  a  hXbXb  mesh  than  with  a  7X7X7 
mesh.  Therefore,  the  lock-in  region  is  larger  with  the  5X5X5  mesh 
(possibly  the  lack  of  any  response  with  the  7X7X7  mesh  indicates 
that  the  disturbance  level  with  a  7X7X7  mesh  faUs  below  the  thresh¬ 
old  amplitude).  We  conclude  that  the  unsteady  results  will  always 
depend  on  the  forcing  amplitude  or  mesh  resolution  in  flows  with  no 
externally  applied  disturbances. 

We’ve  discused  the  oscillation  amplitude  and  frequency,  can  we 
learn  anything  from  the  phase  between  different  probes?  In  the  high- 
density  domain  (see  figure  20),  probes  62,  78,  and  88  oscillate  in  all 
three  directions  while  probes  48  and  120  show  weak  or  no  oscillations. 
In  all  directions,  probes  62  and  78  oscillate  with  approximately  a 
180“  phase  shift,  and  probes  78  and  88  oscillate  with  a  90“  phase 
shift.  The  90“  phase  shift  occurs  at  a  frequency  of  0.054  cycles  per 
time  corresponding  to  a  time  shift  of  4.6  units  between  peaks.  The 
streamwise  distance  between  probes  78  and  88  equals  approximately 
3.0  length  units  in  the  high-density  domain;  therefore,  the  streamwise 
phase  speed  equals 

C’^  =  |^  =  0.65.  (13) 

The  Blasius  velocity  at  the  elevation  of  these  probes  equals  0.62, 
further  evidence  that  the  disturbances  propagate  with  the  local  con¬ 
vection  velocity.  A  disturbance  with  this  frequency  and  phase  speed 
corresponds  to  a  wavelength  of 

A  =  ^  =  12.1.  (14) 

This  wavelength  equals  twice  the  domain  length,  and  the  90“  phase 
shift  between  probes  3.0  length  units  apart  corresponds  to  one  quar¬ 
ter  of  a  wave  with  a  wavelength  of  12.0  length  units. 

How  can  a  disturbance  with  wavelength  equal  to  twice  the  do¬ 
main  length  exist  in  a  domain  with  periodic  boundary  conditions? 
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As  in  all  numerical  simulations,  disturbances  must  satisfy  the  bound¬ 
ary  conditions.  This  may  be  a  restriction  in  one-dimensional  com¬ 
putations;  however,  periodic  boundaries  may  not  restrict  the  dis¬ 
turbance  character  in  two-  or  three-dimensional  domains  when  the 
disturbances  display  only  narrow  spanwise  or  vertical  extent.  When 
a  small  disturbance  travels  downstream  at  an  angle  oblique  with  re¬ 
spect  to  the  direction  of  periodic  boundary  conditions,  or  when  it 
rises  toward  the  free  stream  while  moving  downstream,  there’s  no 
need  for  the  disturbance  to  match  inlet  boundary  conditions  upon 
reaching  the  exit — it  moved  relative  to  its  original  position.  Hence, 
the  disturbance  wavelength  may  be  longer  or  shorter  than  the  do¬ 
main’s  length.  Evidently,  this  is  what  occurs  with  small  roughness 
elements. 

The  oscillations  change  in  the  low-density  domain:  probes  48  and 
120  oscUlate  more  vigorously  than  in  the  high-density  domain,  and 
the  phase  shifts  between  different  probes  no  longer  remain  equal  in 
the  three  directions  (see  figure  21).  In  the  vertical  direction,  the 
phase  shift  between  probes  88  and  78  approximately  equals  5“  (cor¬ 
responding  to  4  time  units);  the  phase  shift  between  probes  62  and 
78  equals  10°  (corresponding  to  8  time  units).  Although  probe  48’s 
results  are  not  shown,  probes  48  and  120  oscillate  nearly  in-phase 
with  a  5°  to  10“  phase  shift  from  the  other  three.  A  zero  phase  shift 
exists  between  probes  62,  78,  and  88  in  the  streamwise  direction.  In 
the  spanwise  direction,  a  small  5“  phase  shift  exists  between  probes 
78  and  88  and  nearly  a  180“  phase  shift  between  these  probes  and 
probe  62.  Clearly  the  simple  structure  evident  in  the  high-density 
case  is  lost  when  the  spheres  are  spaced  further  apart. 

Can  we  classify  the  various  waves  appearing  in  the  different  cases? 
Instabilities  in  periodic  domains  involve  either  standing  or  traveling 
waves.  Since  the  standing  wave’s  phase  speed  equals  zero  (allowing 
independence  of  frequency  and  wavelength)  a  flow  field  dominated  by 
a  single  standing  wave  reveals  signals  from  probes  located  at  different 
points  with  either  a  0“  or  180“  phase  shift.  On  the  other  hand, 
a  traveling  wave  must  must  satisfy  =  /A,  where  the  frequency 
equals  /  and  the  wavelength  equals  A.  In  this  case,  any  phase  shift 
between  probes  is  possible. 

With  the  periodic  pattern  of  spheres,  the  recirculation  zone  be¬ 
hind  each  element  at  a  Reynolds  number  below  that  associated  with 
vortex  shedding  contains  a  standing  wave;  the  entire  separation  zone 
oscillates  at  a  single  frequency.  This  motion  influences  those  particles 
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outside  the  recirculation  zone,  but,  since  they  fall  in  a  region  domi¬ 
nated  by  convection,  the  stamding-wave  description  would  not  apply, 
and  the  flow  contains  a  traveling  wave.  A  probe,  however,  measures 
the  same  frequency  within  the  recirculation  zone  or  outside  in  the 
convection  zone;  although,  the  two  signals  differ  in  phase — phase 
speed  equals  zero  within  the  recirculation  zone,  and  is  finite  outside. 

The  vortex-shedding  mode  at  higher  velocity  contains  no  stand¬ 
ing  waves.  In  this  regime,  shed  vortices  convect  downstream,  and 
numerical  probes  located  at  different  positions  in  the  wake  measure 
a  phase  shift  between  signals  corresponding  to  the  time  for  shed 
structures  to  pass. 

In  the  high-density  domain  the  180^^  phase  shift  between  probes 
62  and  78  indicate  a  standing  wave  in  the  horizontal  direction,  while 
the  90°  shift  between  probes  78  and  88  indicates  a  traveling  wave  in 
the  streamwise  direction. 

In  the  low-density  domain  nearly  all  the  signals  from  the  vari¬ 
ous  probes  show  smadl  phase  shifts  corresponding  to  traveling  waves. 
Only  the  spanwise  signal  from  probes  62  and  88  show  the  character¬ 
istic  180°  phase  shift  associated  with  standing  waves.  The  non-zero 
but  small  phase  shifts  correspond  to  a  rapid  traveling  wave;  the  180° 
shift  may  indicate  a  standing  wave  or  a  traveling  wave  with  wave¬ 
length  equal  to  the  domain  length. 

Apparently,  the  high-density  domain’s  disturbances  behave  like 
standing  waves,  and  the  low-density  domain’s  disturbances  behave 
more  like  traveling  waves.  This  indicates  a  prominence  of  vortex 
shedding  in  the  low-density  domain  and  vortex  waving  in  the  high- 
density  domain. 

5.5  Two  cell  behavior 


What  happens  to  the  characteristic  frequency  after  adding  an¬ 
other  cell  to  the  domain?  Two  possibilities  exist:  the  oscillation 
frequency  drops  by  a  factor  of  two,  as  a  sub-harmonic  oscillation 
with  wavelength  equal  to  twice  the  one-cell  length  occurs;  or  the 
frequency  remains  the  same. 

In  the  first  scenario,  the  periodic  boundary  conditions  must  in¬ 
fluence  the  flow  field,  changing  the  number  of  ceUs  probably  changes 
the  solution,  and  the  individual  roughness  elements  do  not  directly 
produce  disturbances  leading  to  oscillations;  instead,  a  shear  layer 
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based  disturbance  takes  place.  In  the  second  scenario,  it  appears 
that  each  roughness  element  contributes  to  the  disturbance,  the  os¬ 
cillation  seems  directly  linked  to  the  spacing  between  elements,  and 
no  solution  variations  occur  with  changes  in  cell  number. 

Although  the  two-cell  results  are  not  shown,  the  two-cell  high- 
density  domain  yields  St*  =  0.49,  and  the  two-cell  low-density  do¬ 
main  yields  Stk  =  0.24.  These  agree  with  the  one-cell  results,  and 
it  appears  that  the  single  cell  adequately  accounts  for  the  important 
physics.  Furthermore,  the  various  phase  shifts  between  signals  from 
different  probes  match  in  the  one-  and  two-cell  cases. 

5.6  Comparison  with  experiments 

Experiments  [Tadjfar  &  Others,  1985]  at  Rck  =  310  and  k/6*  = 
1.6,  w'*h  an  array  of  spherical  roughness  elements  similar  to  the 
standard  density  case  used  in  this  study,  reveal  Stk  =  fk/Uk  =  0.08 
using  the  undisturbed  Blasius  velocity  and  0.23  using  the  measured 
velocity  at  the  top  of  a  sphere.  These  results  compare  favorably 
with  Stk  =  0.23  with  isolated  spheres  [Mochizuki,  1961],  but  dif¬ 
fer  from  the  Stk  =  0.33  numerical  results.  This  disparity  cannot 
be  resolved  without  additional  experimental  results  using  a  different 
roughness  density.  The  numerical  and  experimental  results  agree  in 
one  area,  however.  The  peak  amplitude  occurred  near  1.5fc  from  the 
waU  in  both  cases,  although  Tadjfar  took  no  measurements  within 
the  roughness  elements. 

6  Conclusions 

Since  identical  oscillation  frequencies  were  obtained  with  peri¬ 
odic  boundary  conditions  applied  across  both  one-  and  two-cell  do¬ 
mains,  the  domain  length  does  not  govern  the  oscillation  frequency; 
instead,  the  primary  instability  associated  with  a  periodic  array  of 
spheres  protruding  from  a  flat  surface  occurs  at  a  frequency  which 
depends  on  the  streamwise  spacing  between  elements  and  the  lo¬ 
cal  velocity  at  the  elevation  of  the  sphere’s  crest.  This  also  means 
that  the  characteristic  frequency  remains  independent  of  the  isolated 
sphere  shedding  frequency,  and  any  similarities  between  the  two  are 
merely  coincidental. 
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This  oscillation  frequency,  called  the  passing  frequency  for  the 
time  required  for  a  fluid  particle  to  pass  from  one  sphere  to  the  next, 
may  also  be  thought  of  as  a  forcing  frequency  due  to  disturbances 
generated  in  the  wake  region  behind  the  individual  spheres.  The 
response  where  the  oscillation  frequency  equals  the  forcing  frequency 
is  called  a  lock-in.  In  the  lock-in  mode,  the  frequency  of  oscillation 
must  lie  near  the  shedding  frequency  of  a  similar  body  located  in  an 
infinite  domain — in  this  case,  the  similar  body  equals  a  single  sphere 
attached  to  a  flat  surface. 

At  Rck  =  175,  no  streamwise  vortices  based  on  the  flow’s  inertia 
exist;  instead,  the  flow  field  contains  regions  of  streamwise  vorticity 
explained  entirely  by  continuity  considerations. 
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13.  ABSTRACT  (Maximum  200  words) 

Results  from  a  direct  numerical  simulation  of  laminar  flow  over  a  flat  surface  with  spherical  roughness  elements  using  a  spectral- 
element  method  are  given.  The  numerical  simulation  approximates  roughness  as  a  cellular  pattern  of  identical  spheres  protruding  from 
a  smooth  wall.  Periodic  boundary  conditions  on  the  domain’s  horizontal  faces  simulate  an  infinite  array  of  roughness  elements 
extending  in  the  streamwise  and  spanwise  directions,  implies  the  parallel-flow  assumption,  and  results  in  a  closed  domain.  A  body 
force,  designed  to  yield  the  horizontal  Blasius  velocity  in  the  absence  of  roughness,  sustains  the  flow.  Instabilities  above  a  critical 
Reynolds  number  reveal  negligible  oscillations  in  the  recirculation  regions  behind  each  sphere  and  in  the  free  stream,  high-amplitude 
oscillations  in  the  layer  directly  above  the  spheres,  and  a  mean  profile  with  an  inflection  point  near  the  sphere's  crest.  The  inflection 
point  yields  an  unstable  layer  above  the  roughness  (where  Uly)  <  0)  and  a  stable  region  within  the  roughness  (where  Ulyl  >  0). 
Evidently,  the  instability  begins  when  the  low-momenfum  or  wake  region  behind  an  element,  being  the  region  most  affected  by 
disturbances  (purely  numerical  in  this  case),  goes  unstable  and  moves.  In  incompressible  flow  with  periodic  boundaries,  this  motion 
sends  disturbances  to  all  regions  of  the  domain.  In  the  unstable  layer  just  above  the  inflection  point,  the  disturbances  grow  w  hile  being 
carried  downstream  with  a  propagation  speed  equal  to  the  local  mean  velocity;  they  do  not  grow  amid  the  low  energy  region  near  the 
roughness  patch.  The  most  amplified  disturbance  eventually  arrives  at  the  next  roughness  element  downstream,  perturbing  its  wake 
and  inducing  a  global  response  at  a  frequency  governed  by  the  streamwise  spacing  between  spheres  and  the  mean  velocity  of  the  most 
amplified  layer. 
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